Hyperspherical explicitly correlated Gaussian approach for few-body systems with 

finite angular momentum 



D. Rakshiti and D. Blumei-^ 

^Department of Physics and Astronomy, Washington State University, Pullman, Washington 99164-2814, USA 

^ITAMP, Harvard-Smithsonian Center for Astrophysics, 
60 Garden Street, Cambridge, Massachusetts 02138, USA 
(Dated: December 27, 2012) 

Within the hyperspherical Iramework, the solution of the time-independent Schrodinger equation 
for a n-particle system is divided into two steps, the solution of a Schrodinger like equation in the 
hyperangular degrees of freedom and the solution of a set of coupled Schrodinger hke hyperradial 
equations. The solutions to the former provide effective potentials and coupling matrix elements that 
enter into the latter set of equations. This paper develops a theoretical framework to determine the 
effective potentials, as well as the associated coupling matrix elements, for few-body systems with 
finite angular momentum L = 1 and negative and positive parity D. The hyperangular channel 
functions are expanded in terms of explicitly correlated Gaussian basis functions and relatively 
compact expressions for the matrix elements are derived. The developed formalism is applicable to 
any n; however, for n > 6, the computational demands are likely beyond present-day computational 
capabilities. A number of calculations relevant to cold atom physics are presented, demonstrating 
that the developed approach provides a computationally efhcient means to solving four-body bound 
and scattering problems with finite angular momentum on powerful desktop computers. Details 
regarding the implementation are discussed. 

PACS numbers: 



I. INTRODUCTION 

Few-body phenomena play important roles across all 
disciplines of physics, including atomic and molecular 
physics, chemical physics, nuclear and particle physics, 
and condensed matter physics. Progress in solving and 
understanding the quantum mechanical fevir-body prob- 
lem has been driven, roughly speaking, by one or more of 
the following three aspects: (i) Using well-established al- 
gorithms, the steady increase of computational resources 
has made it possible to tackle problems that were impos- 
sible to tackle a decade or even just a few years ago. (ii) k 
number of model systems have been investigated analyt- 
ically, semi-analytically or numerically, providing crucial 
insights into some of the low-energy processes that gov- 
ern the few-body dynamics, (in) More efficient numerical 
schemes that are not only applicable to the three-body 
problem but also to four- and higher-body problems have 
been developed. 

This paper extends the correlated Gaussian hyper- 
spherical (CGHS) or hyperspherical explicitly correlated 
Gaussian (HECG) approach [J-il. In earlier work, von 
Stecher and Greene [l-Q considered three- and four-body 
systems with vanishing angular momentum L and posi- 
tive parity 11. Here, we extend the approach to systems 
with finite angular momentum L. Although the overall 
scheme developed for systems with = 0^ symmetry 
carries over to systems with finite angular momentum, 
the determination of compact expressions for the matrix 
elements associated with finite angular momentum states 
is significantly more involved than that for states with 
= 0+ symmetry. 

The HECG approach provides an efficient numerical 



scheme for solving few-body problems. It is a basis set 
expansion type approach, which combines elements of the 
aspects (i)-(iii) mentioned above. In particular, the use 
of hyperspherical coordinates [sl-fTsj within the frame- 
work of explicitly correlated Gaussian (EGG) basis func- 
tions [l6| allows us to take advantage of the machinery 
developed for bound state calculations while at the same 
time enables us to describe the scattering continuum. 
Although significant progress has been made [l|, |^, [TtI - 
[20| . in general, the determination of scattering quantities 
is significantly more involved than that of bound state 
quantities, and the four-, five- and higher-body scatter- 
ing continua are comparatively poorly understood. Thus, 
the framework developed in this work for systems with 
1~ and 1"*" symmetry provides a promising step forward. 

The HEGG approach is quite general and applicable 
to a wide range of few-body systems. As an applica- 
tion, we consider the four-particle system consisting of 
three identical fermions and an impurity whose mass is 
lighter than that of the majority species. We assume in- 
terspecies short-range s-wave interactions and investigate 
the system properties of the energetically lowest-lying 1+ 
state as a function of the mass ratio k between the ma- 
jority particles and the impurity particle. These finite 
angular momentum states are interesting since univer- 
sal four-body bound states have been predicted to exist 
if the two-body s-wave scattering length is positive and 
K > 9.5 dJ. Moreover, for 13.38 < k < 13.61, the (3, 1) 
system with 1+ symmetry has been predicted to support 
four-body Efimov states in this mass ratio regime, 
three-body Efimov states are absent ^23.',-26j. This work 
determines and interprets the hyperangular eigen value 
of the (3, 1) system with infinitely large interspecies s- 
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wave scattering length in the hmit that the hyperradius 
is much larger than the range of the underlying two-body 
potential. 

The remainder of this paper is organized as follows. 
Section introduces the system Hamiltonian and the 
hyperspherical framework, while Sec. IIIII introduces the 
ECG basis functions used to expand the hyperangular 
channel functions. Section IIVI discusses the matrix ele- 
ments, applicable to any n, needed to calculate the effec- 
tive hyperradial potential curves and associated coupling 
matrix elements. Details regarding the numerical imple- 
mentation of the HECG approach and a set of proof- 
of-principle calculations for the four- and five-particle 
system are discussed in Sec. |Vl Section |Vl] applies the 
HECG framework to the (3, 1) system with 1+ symme- 
try and diverging interspecies s-wave scattering length 
for various mass ratios k. Lastly, Sec. IVIII summarizes 
and concludes. Details regarding the derivation of and 
final results for the fixed hyperradius matrix elements 
are presented in three appendices. Appendix 1X1 defines a 
number of auxiliary quantities that depend on the sym- 
metry considered. Appendix IB] outlines exemplarily how 
to derive the matrix elements for the three-body system 
with 1~ symmetry. Appendix [Cl summarizes our expres- 
sions for a number of quantities that enter into the final 
equations for the matrix elements; these equations apply 
to all symmetries considered in this paper. 



II. SYSTEM HAMILTONIAN AND 
HYPERSPHERICAL FRAMEWORK 
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where fj,j denotes the mass associated with the j*'^ Jacobi 
vector. 



rrik rrij+i 



for j = 1, • • • ,n - 1 



and 



Mn = Wfc. 
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The generalization to n > 5 is straightforward. By def- 
inition, the n"^ Jacobi vector coincides with the "mass- 
scaled" center of mass vector of the n-particle system. 
Although the mass-scaling is not needed to separate ofi^ 
the center of mass motion, the use of mass-scaled Jacobi 
vectors — as opposed to the use of non-mass-scaled Jacobi 
vectors — simplifies the derivation of the fixed- i? matrix 
elements for the ECG basis functions (here, R denotes the 
hyperradius; see below and Sec. lIV[) . The Hamiltonian H 
can now be written as a sum of the relative Hamiltonian 
H^ci and the center of mass Hamiltonian Hem, 



We consider an n-particle system with position vectors 



fj described by the Hamiltonian H, 



(1) 



where mj denotes the mass of the j particle. The in- 
teraction potential Vint is written as a sum of two-body 
potentials Vjk{rjk), 



Vint = 



(2) 



where fjk — Vj —rk {fjk — l^jfcD- To separate off the cen- 
ter of mass degrees of freedom, we define n mass-scaled 
Jacobi vectors pj , 



fc=i 



(3) 



The elements Tjk form anxn matrix. The explicit forms 
for n — 5 and n — A read 
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with 



H = H,, 



fc2 

2 ^ Pj 



and 



ticm - — 2 
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In the following, we seek solutions to the relative 
Schrodinger equation 



(12) 



i.e., we seek to determine V'e Sind E. The energy E can 
be negative or positive, i.e., we consider both bound state 
and scattering solutions. 

We employ the hyperspherical coordinate approach [1, 
0, H, [iB|, which has proven to provide critical physical 
insights that, in some cases, are more difficult or even 
impossible to unravel in alternative approaches. The so- 
lution to the relative Hamiltonian is divided into two 
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steps: (i) the solution of a Schrodinger like equation 
in the hyperangular coordinates and (ii) the solution of 
a Schrodinger like equation in the hyperradial coordi- 
nate. More specifically, the idea is to expand the relative 
wave function V'_b(pi, • • • , Pn-i) in terms of a complete 
set of hyperangular channel functions 57) that de- 

pend parametrically on the hyperradius R and hyperra- 
dial weight functions FyE{R) [1,0, Sill, 



ill 



R 



-(3n-4)/2 



Here, R denotes the hyperradius. 



n-l 

R' = Y.pI 



(13) 
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which has, as the components of pj, units of "mass^/^ 
times length". The mass-scaled hyperradius R can be 
related to the "conventional unsealed hyperradius" by 
pulling out a factor of ^/Jl, where /x is the hyperradial 
mass. In Eq. (IT3)) . 51 collectively denotes the 3ri — 4 hy- 
perangles. The hyperangles fl can be defined in differ- 
ent ways (see Sec. IIVI for the definition employed in this 
work) . 

The channel hmctions $^(i?;il) form a complete set 
in the (3n — 4)-dimensional Hilbert space associated with 
the hyperangular degrees of freedom 3, 7, 8, 15], 



[$,.(i?; r!)]*$,(i?;r!)d-^"-4ll = 5,^, 



(15) 



The <i>,y(i?; are chosen to solve the fixed- i? hyperangu- 
lar Schrodinger equation 

ifadia + Vn,t{R, ^)\ $.(i?; Vt) - U,{R)^,{R- 9), (16) 

where 

i?adia=ro+Kfi(i?) (17) 



with 
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The grandangular momentum operator A [l,[l3] accounts 
for the kinetic energy associated with the hyperangular 
degrees of freedom. For our purposes, it proves advanta- 
geous to define Tn through 



where 



2 R^'^-^dR OR' 
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Inserting Eq. (IT3|) into Eq. (IT2l) . we find that the weight 
functions F^e (R) and the relative energy E are obtained 
by solving a set of coupled hyperradial equations 
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where the coupling term Vc^u is given by 

VcAR) = 

^ r { r>\ 
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with 



P.AR) = ^ ImR-Ar^-^^^^^d'-'^ (24) 
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To reiterate, the hyperspherical framework consists of 
two steps: In the first step, the (3n — 4)-dimensional 
hyperangular Schrodinger equation is solved, yielding 
Ui,{R), Puu'iR) and Qyv'{R)- In the second step, the 
coupled set of one-dimensional hyperradial equations is 
solved, yielding F^e{R) and E. This paper focuses pri- 
marily on solving the hyperangular Schrodinger equation. 
Expressions for the relevant matrix elements, valid for 
any n, are derived and applications to systems with n = 4 
are presented. 



III. FUNCTIONAL FORM OF THE BASIS 
FUNCTIONS 

To solve the hyperangular Schrodinger equation, we 
expand the channel functions $!y(i?;f2) for fixed R in 
terms of ECG basis functions ik{A^^\u^^\i62\x)\R^ 



Xj\R, 



(26) 



k=l 



where x collectively denotes the 3n — 3 Jacobi vectors, 
X = (pi, • • • ,/3„_i). In Eq. ([26|l . the notation "|i?" indi- 
cates that Tpk is evaluated at a fixed hyperradius R. The 
(n — 1) X (n — 1) dimensional matrices -A^'^-' are symmet- 
ric and positive definite. The n{n — l)/2 independent 
elements of -A*-'^-' are treated as variational parameters 
of the /c'^ basis function. The elements of the (n — 1)- 
dimensional vectors u'^i^ and u':2^ are also treated as vari- 
ational parameters. The optimization scheme employed 
to determine the values of these variational parameters 
is discussed in Sec. |Vl In Eq. I^E^, S denotes an opera- 
tor that imposes the proper symmetry under exchange of 
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identical particles. For the three-body system consisting 
of two identical fermions and an impurity, e.g., we have 
S — 1 — Pi2- For the four-body system consisting of three 
identical fermions and an impurity (see Sec. lVip . we have 

S — 1 — Pi2 — Pl3 — P23 + ^"12^23 + -Pl3-P32- 

The functional form of the basis functions depends on 
the symmetry considered. We consider ECG basis 
functions of the form 16, 27-31] 



1p{A,Ui,U2,x) = Nlh X 



vi\''\v2r[Yi,{vi)(g,Yi,{v2)]LM^ exp 



(27) 



which can be used to describe states with 
0+, 1+, 1^, 2+, • • • symmetry but not states with 
0~ symmetry (32j . The three-dimensional vectors vj, 
j — 1 and 2, are defined through 



n-l 
k=l 



(28) 



where Uj^k denotes the k component of the vector Uj. 
We denote the elements of the vector Vj by Wj^i, Vj^2 and 
Vj^3. In Eq. (P7)) . the notation \Yi^ ® Yi^]lml indicates 
that the two spherical harmonics Yi-m^ are coupled to a 
function with total angular momentum L and projection 
quantum number M^, and N^yi denotes a normalization 
constant. 

In the following, we write the basis functions out for 
Ml — 0; in writing the basis functions for a specific 
symmetry, we choose the normalization constant A^lii hi 
Eq. ([77)1 conveniently. Throughout this paper, we restrict 
ourselves to states with 0"'",1~ and 1+ symmetry. The 
= 0+ basis functions, obtained by setting li and I2 to 
0, are independent of ui and U2 (or, equivalently, of vi 
and V2), 



ipiAi x) — exp 



(29) 



The = 1 basis functions, obtained by setting li to 1 
and Z2 to 0, depend on ui but not on U2, 



ui, x) = 3^/^1)1,3 exp 



'Ax 



(30) 



Lastly, the = 1+ basis functions, obtained by setting 
li and ^2 to 1, depend on both ui and M2, 

■4'{A,ui,U2,x) = 



21/2 



(Wl,2«2, 



I'l, 1^2,2 J exp 



X Ax 



(31) 



The next section presents relatively compact expres- 
sions for the fixed- i? matrix elements between two basis 
functions ip = ip(A,ui,U2, x) and ip' = ip(A ,Ui,U2,x). 
Throughout, we assume that tp and ip' are characterized 
by the same L, Ml and 11 quantum numbers. For sys- 
tems with finite angular momentum, neither the fixed- i? 
overlap matrix element nor the fixed- i? matrix elements 
for Tq , Pjyi,/ , Q 1,1,1 , and Vint have, to the best of our knowl- 
edge, been reported in the literature. 



IV. MATRIX ELEMENTS FOR 
HYPERSPHERICAL EXPLICITLY 
CORRELATED GAUSSIANS 

We introduce the short-hand notation 

(V''IV')U = 

mA\ u'l , "2 , a?) U] > (A, «i , W2 , x) I Rd^^'-^n, (32) 



— / [V'(A',u'i,U^,f)|fl,] 



- /[v(a',u;,4,^)Ih]^ 

and 



.,dtp{A, ui,u2,x) 



dR 



73n-4 



17,(33) 



(^-'1^1^)1^ 



■d'^lpiA, Ul,U2,x) 



dR2 



^ / [^(^',w'i,4,^)|fl]*7^oV'(A"l,W2,f)Ud3"-4d- 



[V'(A ui,u2,x)\R]*Tn^{A', u[,u'2, f)Ud=^"-4ri(35) 



As in Refs. [1H3|, Eq. ([55)) explicitly symmetrizes the ma- 
trix element associated with the hyperangular kinetic en- 
ergy. 

To evaluate the fixed-i? matrix elements defined in 
Eqs. we introduce a new set of coordinates y, 

2/ = (yi,--- ,2M-i), through y = {HbYx. The matrix 
is chosen such that the matrix (U p)^ BU where 



B = A + A', 



(36) 



is diagonal with diagonal elements • • • , Pn~i, i.e., such 
that x^Bx — f^iVj- follows that the arguments 

of the exponentials in the integrals defined in Eqs. (j32p - 
(1351) reduce to quadratic forms. Since the coordinate 
transformation from x to y is orthogonal, we have (i) 
i?2 ^ ^Jji yfi and (it) / • • • = / . . . d3("-i)y. 

To perform the integration over fl, we need to specify 
the 3n — 4 hyperangles. Following Refs. [iHl, [13, HB], we 
define 2(n — 1) angles as the polar and azimuzal angles 
"dj and ipj of the n — 1 vectors yj. The remaining n — 
2 hyperangles 71, • • • ,7ti-2 are defined in terms of the 
direction of the (n — l)-dimensional vector s, where s = 
{yi, ■ ■ ■ ,yn-i) and j/j = \yj\. Specifically, we write 

yi = i?sin7isin72 sin7„_2, 

j/2 = i?cos7isin72 sin7„„2, 

y^ = R cos 72 sin 73 sin 7„_2 , 

* * ' 7 

yn-2 = i?cos7„_3sin7„_2, 

yn-i = i?cos7„_2, (37) 
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where 7^ € [0,7r/2]. This restriction on the range of the 
angles ensures that the yj are non-negative. Correspond- 
ingly, we have [s^, l35j | 



■■■d 



'3ri-4 



n = 



n-2 




7fc cos^ 7fcd7fc , (38) 



fc=i 



where (Pijj denotes the "usual" angular piece of 
the volume element in spherical coordinates, cPyj = 
sin -dj ddj dipj . 

In general, the expressions for the matrix elements de- 
fined in Eqs. ([5^ - (|55|) depend on the functional form of 
the basis functions ^p{A,ui,U2,x) (see Sec. IIII|) . For the 
basis functions defined in Eqs. (|29)) - (|3T|) . the integrations 
involving the angles "dj and ipj (j = 1, • • • , n — 1) can be 
performed analytically, yielding 



tional forms. We write 



+ ^ [, 

n-l 



22) 2 2 

VjVk 



,(44) 4,4 



y.yk 



k>j=l 
n-l 



V 4 ,(26) 2 6 

j=l,k=l,k^j 



,(222) 2 2 2 

d-),k,iyjVkyi 



,(224) 2,,2„4 



j^k/yjykyi 

{2222) ,,2„2„2,,2 



E 

A;>j=i;i7^i,fc 

E (43) 

k>j—l;l>j;'m>l\m^l^k^j 



ii^'Mn = (47r)"-i / /<,(s)exp ( -^^fty,^ 



n 



sm 



3fc-l 



7fcC0S^7fcrf7fe , (39) 



(^'|P|^)|« = - 



;i2(47r)^ 



ft2(47r)"- 



and 



(^-'1^01^)^ = - 



;i2(47r) 



/p(.s)exp 



/Q(s)exp 



In writing /o(s), we dropped terms that contain odd pow- 
ers of yj since these terms average to zero when inte- 
grating over the remaining hyperangles. The quantities 
/p(s), /q(s) and fn{s) are obtained by replacing the d's 
in Eq. p3|) by p's, q's and 6's, respectively. The super- 
and subscripts of the d-, p-, q- and 6-coefficients indicate 
the polynomial in the y's that the coefficients are asso- 
ciated with. The d-, p-, q- and 6-coefFicients depend on 
the symmetry of the wave function, and are listed in Ap- 
pendix]^ for states with = O"*", 1^ and 1+ symmetry. 
It should be noted that, depending on the symmetry and 
number of particles, a varying number of the d-, p-, q- and 

) 6-coefficients vanish. Appendix |B] exemplarily illustrates 
(,40) how to obtain Eq. dSH]) for the three-body system with 
= 1^ symmetry. We emphasize that Eqs. 
together with the expressions given in Appendix |X1 apply 
to any number of particles. For states with L > 1 and 
= 0~, the definitions of the d-, p-, q- and 6-coefficients 
contained in /o, /p, /q and /n change and polynomials in 
the y's of higher power than those considered in Eq. p3)) 
may appear. 



i=i 



sm 



3fc- 



^7feC0S^7fcd7fc (,41) 



^ n — 1 

/o(s)exp I - 2E^jy^ 



The integration over 71 in Eqs. can also 

be done analytically. To perform this integration, we 
recognize that the hyperangle 71 enters into yi and 
y2 but not into yj with j > 3. Using Eq. (|37p . 
we write yj = i?^ sin^ 71 i7(72, 7,1-2) and y| = 
R'^ cos'^ jiH{'j2, - ■ ■ ,7n-2), where H{j2,--- ,7^-2) = 1 
for n = 3 and 



sm 



3k- 



SfcCOS^7fcd7fe (42) 



Hh2 



, 7n-^2) = sm 72 X • • • X sm 7„ 



(44) 



The matrix elements {^'IPMr, {^'IQMb^ 

and {■ilj'\Tn\ilj)\R, Eqs. (gUl-dlll), have been written such 
that /o(s), /p(s), /q(s) and fn{s) have analogous func- 



for n > 3. In the following, we suppress the dependence 
of H on the angles 7^ {j > 2) and rewrite /o(s) such that 
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the dependence on 71 is "isolated" , 

fo{s) = SCOQ + SC2qHR^ SVC? 71 + SC02HR^ cos^ 71 + 

sc4o{HR'^f sin^ 71 + sco4{HR^f cos^ 71 + 
sceoiHR^f sin^ 71 + scoe{HR^f cos^ 71 + 
SC22{HR^Y siii^ 7i cos^ 71 + 
sCi4{HR^Y sin'' 71 cos* 71 + 
sc24(-ffi?^)^ sin^ 7i cos'' 71 + 
sci2{HR^)^ sin* 71 cos^ 71 + 
sc2q{HR?Y sin^ 71 cos^ 71 + 
sc62(-ffi?^)*sinSicos2 7i- (45) 

The coefficients scjk depend on the liyperangles 7; with 
I > 2 and the d-coefficients, and are defined in Ap- 
pendix [C] The subscripts j and k of the sc-coefficients 
denote respectively the powers of sin 71 and cos 71 that 
the coefficients scjk are associated with. Using Eq. (1451) . 
we find 



tt/2 



/o(s)exp 



1 



sin^ 71 cos^ 71^71 = 



Ml 



MO 

c 



/2(C) 



c 



(46) 



where /i and I2 denote Bessel functions and C is defined 
through 



c 



^{P,-f32)HR^ 



(47) 



The quantities Mi and M2 can be written in terms of 
RH^, C and the sc-coefficients; Mi and M2 depend on 
727 ■ ■ • ,7n-2 since H and the sc-coefficients depend on 
these angles. Explicit expressions for Mi and M2 are 
given in Appendix [Cl Using Eq. (|46)) in Eq. ([39)) . the ma- 
trix element {ip\ip)\R is known fully analytically for n = 3 
and reduces to a (n — 3)-dimensional integral for n > 3. 
For n — 4 and 5, the remaining one- and two-dimensional 
integrations can, as discussed in Sec.|Vl be performed nu- 
merically with high accuracy. Expressions P5|) and pS)) 
also apply to /p(s), /q(s) and /o(s) if the ^-coefficients 
are replaced by the p-, q- and 6-coefficients, respectively. 

For = 0+, our results obtained using the above 
expressions agree with those reported in Refs. P, [2| for 
n = 3 and 4. Motivated by our desire to express the 
various matrix elements for different number of particles 
n and symmetries in a unified framework, we adopted 
a notation that differs notably from the notation adopted 
in Refs. [H-i,!!^. 

We now turn to the evaluation of the interaction matrix 
element. We define 



/ 



u'i,u'2,x)\R]*Vki{fkimA, ui,U2, a;)|Krf3"-4q;.48) 



If the two-body potential Vki{rki) is parameterized by a 
spherically symmetric short-range Gaussian Vg{rki) with 
depth Dki and range ro,fc/, 



-Dki exp 



V2\ 



ro,ki 



(49) 



then {i;'\Vg{rki)\ip)\R is equivalent to -Dki{ip'\'fp)\R. if 
A and A' are replaced hy A + iii^'''V(2''o,fci) and A' + 
W^''''^ /{2rl f^i), respectively. Here, the matrix w}-^'^^ is 
defined as 



(50) 



where the vector cJ^*^') provides the transformation from 
the Jacobi coordinates x to the interparticle distance vec- 
tor rki, 



ru = 



(51) 



(52) 



It follows that we can use Eq. ([55]) [see also Eqs. P5)) - 
(H71) ] if C is replaced by C}-^^^ and if the /3j are replaced 
by /Sj*^'-' . Similar to C, C^*^''' is defined through 

CC^O ^i(/3f')-/3f))iJi?2 
and the /jj'^'' denote the eigenvalues of the matrix B_ 



V. IMPLEMENTATION DETAILS OF THE 
HECG APPROACH AND 
PROOF-OF-PRINCIPLE APPLICATIONS 

As discussed in the previous sections, the determi- 
nation of the effective hyperradial potential curves and 
coupling matrix elements requires that the hyperangu- 
lar Schrodinger equation be solved for several hyperradii 
R. For each fixed i?, the determination of the linear 
and non-linear variational parameters is accomplished 
following approaches very similar to those employed in 
the "standard" (non fixed R) EGG approach For 
a given set of basis functions, and thus for a given set 
of non-linear variational parameters, the expansion coef- 
ficients Cfe, k = 1, • • • , A^fc [see Eq. p6|) ]. are determined 
by diagonalizing the generalized eigenvalue problem de- 
fined by the fixed- i? Hamiltonian matrix and the fixed- i? 
overlap matrix. According to the generalized Ritz vari- 
ational principle, the Nf, eigenvalues provide variational 
upper bounds to the exact eigenvalues of the hyperangu- 
lar Schrodinger equation. 

The non-linear variational parameters are collected in 
A^^\ u^i^ and Mj'^^ where k = I,-- - ,A^fc, and deter- 
mined using the stochastic variational approach (STj . i.e., 
through a trial and error procedure. For concreteness. 
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we consider the situation where we aim to determine the 
energeticahy lowest lying hyperangular eigenvalue Uq (R) 
for a given R value. We start with a small basis set 
(typically consisting of just one basis function) . To add 
the next basis function, we semi-randomly generate 
trial basis functions (A't is typically of the order of a few 
thousand), yielding trial basis sets. We determine 
the eigenvalue for each of these trial basis sets and choose 
the trial basis set that yields the smallest eigenvalue as 
the new basis set. Following the same selection process, 
we continue to enlarge the basis set one basis function 
at a time. This procedure is repeated till the basis set 
is sufficiently complete and the desired accuracy of the 
energetically lowest lying eigenvalue is reached. The opti- 
mization of excited states proceeds analogously. If nearly 
degenerate states exist, it is advantageous to simultane- 
ously minimize the eigenvalues of multiple states. 

Since the trial and error procedure "throws out" most 
of the trial basis functions generated, the resulting basis 
set is typically comparatively small. For the four-body 
systems discussed below, e.g., we achieve convergence for 
iVb of the order of 100. Moreover, we have found that a 
carefully constructed basis set avoids a number of prob- 
lems that can arise from the fact that the basis func- 
tions are not orthogonal. In particular, the minimiza- 
tion scheme that underlies the trial and error procedure 
tends to select basis functions that have relatively small 
overlaps among each other, i.e., that are "fairly linearly 
independent" . In some cases, however, the trial and er- 
ror procedure does not fully eliminate numerical issues 
arising from the linear dependence of the basis functions. 
Thus, we add another check and reject a given trial basis 
function V'T if its overlap with one or more of the basis 
functions already selected is too large, i.e., if the quantity 
('0t|V'/c)|/j is larger than a preset value e, where we nor- 
malize ipT and tpk such that (i/'T|V'T)|i?, — (V'fclV'fc)!/? = 1- 
We have used e = 0.9 or 0.95 in most of the calculations 
reported below. A similar criterion is employed in the 
HECG approach of Ref. Q and in the standard EGG 
approach |16|. 

We now discuss the determination of the non-linear 
parameters that characterize the basis functions. The 
selection of the parameters of the matrices A^'"'^ is guided 
by physical considerations. The fact that the basis func- 
tions have to be "compatible" with the value of the hy- 
perradius considered implies that the parameters of the 
matrices have to be chosen such that the quantity 
J2'j<ii^''ji ^)'^ is of the order of R^ / ^. The width parame- 
ters are related to the parameter matrix A^*'-' via 



E 

3<l 



(53) 



Moreover, the basis functions have to govern the dynam- 
ics that occurs at the length scale of the two-body inter- 
action potential. Gorrespondingly, we consider different 
types of basis functions: The first type is characterized 



by all cf^i^ being of the order of R/{n^JJi)] the second 

type is characterized by one cf^j^ being of the order of 
the range of the underlying two-body potential (we use 
To to denote the smallest of the J^o.ji's) and all other d -i 
being of the order of R/{n^JJl)] the third type is charac- 

(k) 

terized by two d -i being of the order of j-q and all other 
d'^^l' being of the order of R/{ny/Ji); and so on. The ex- 
act values of the d^-J^ are chosen stochastically from pre- 
defined parameter windows, which are chosen according 
to the above considerations. We have found that the 
convergence of the hyperangular eigenvalues for strongly 
interacting systems with large R/{y/Jlro) depends quite 
sensitively on the choice of the parameter windows. As 
in Ref. 1] , we allow for basis functions with positive and 
negative widths. The elements of the parameter vectors 
M^'^^ and u^2^^ are chosen to lie between —1 and 1. 

The quantity |^|, see Eq. (|47)) . takes on large values 
if |/3i| > 1/32 1 or 1/32 1 > This situation arises 

quite frequently if the hyperradius R is much larger than 
^/JItq and causes numerical difficulties when evaluating 
the Bessel functions. We mitigate these difficulties as fol- 
lows. For concreteness, we consider the n = 4 case and 
the integral involving /i(^). For large |^|, we write the 
integral over the hyperangle 72 [see Eqs. (l39l) and (|46p ] 
as 



/'r'exp 
Jo 



C ^ cxp 



-i(/3i + /32)i?2(l-x2)-i/33i?V 

M,{x)^-^il~x^fx^dx^ 
(27r)-i/2 , 

-i min(/3i, /32)i?'(l - x^) - )^^zR^x' 



15 



105 



|^|3/2 8|^|5/2 128|^|V2 1024|C|9/^ 

M^(x){\~x^fx^dx{h4) 



where x — cos 72 [implying -^(72) = 1 — x^]. We use the 
right hand side of Eq. ([M)) when |^| > 50. An analo- 
gous expansion is done for the part of the integrand that 
involves /2(0- 

Although M\ depends on the overall behavior of 
the integrand in Eq. (j54p is determined by the exponen- 
tial. In particular, depending on the signs of min(/3i, /32) 
and /33, the integrand can be sharply peaked at a; w 
or X ~ 1. To perform the integration over x (i.e., the 
integration over the hyperangle 72) numerically, we di- 
vide the integral into three sectors: the left, middle and 
right sectors. The sector boundaries Xmid,! and Xniid,2 
are determined dynamically depending on the behavior 
of the integrand. The default values are a;mid.i — 0-3 and 
2:mid,2 = 0.7. Howcvcr, if /Ssi?^ > b or min(/?i, /32)i?^ < 
— (5, we choose Xinid.i = l/\/max(/33_R2^ | min(/3i, /32)|i?^). 
Similarly, if fi^R? < —S or min(/3i, /32)-R^ > 5, we choose 
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a;mid,2 = 1 - l/v/niax(|^3|i?2, miii(/3i, /32)i?^). Our cal- 
culations reported below use 6 = 20. The numerical in- 
tegration for each sector is performed using a standard 
Gauss-Laguerre integration scheme Typically, we 

choose the same order A'ordor for all three sectors [the 
value of Nb depends on the ratio R/ {^/Jlro)]. We note 
that the integration scheme employed has not been opti- 
mized carefully and can possibly be refined further. For 
the five-body system, appropriate generalizations are in- 
troduced. 

To demonstrate that the developed framework works, 
we consider four-body systems consisting of two identical 
spin-up fermions and two identical spin-down fermions. 
Such systems can be realized experimentally by occu- 
pying two different hyperfine states of ultracold atomic 
^Li or samples |39l - l4l| . We assume that the identical 
fermions do not interact. This is a good assumption since 
s-wave interactions between identical fermions are forbid- 
den by the Pauli exclusion principle and p-wave interac- 
tions are, in most experimental realizations, highly sup- 
pressed by the threshold law. Furthermore, we assume 
that the interspecies interactions have been tuned so that 
the interspecies free-space s-wave scattering length as is 
infinitely large. This regime is referred to as the unitary 
regime and can be realized experimentally by applying an 
external magnetic field in the vicinity of a Fano-Feshbach 
resonance [42|. In our calculations, we describe the inter- 
species two-body interactions using a Gaussian two-body 
potential with depth D and range tq adjusted such that 
the two-body potential supports a single zero-energy s- 
wave bound state. We calculate the energetically lowest- 
lying hyperangular eigenvalue Uq {R) for various ^/Jlro / R 
values. 

To present our results, we rewrite Uo{R) in terms of 
the quantity so{R), 



UoiR) 



hmsoiR)]^~l/A} 
2i?2 



(55) 



The scaling introduced in Eq. (1551) is motivated by the 
fact that so(i?) becomes independent of R in the non- 
interacting limit and if ro <C \as\ [ij. The latter regime 
is realized if the s-wave scattering length diverges and if 
the quantity y/JIro is much smaller than R. 

Circles in Fig. [T] show the quantity so.unit(-R) for the 
(2,2) system at unitarity as a function of ^JJir^/R for 
(a) — O'^ symmetry, (b) = 1~ symmetry, and (c) 
= 1+ symmetry. The basis set extrapolation error is 
smaller than the symbol size. Dotted lines show a fit to 



the data using the fitting function Sq 



ZR 



^//Iro /R. The coefficient Sq ^^j^ is found to be 



where x 

2.510(1), 4.600(3), and 4.081(3) for = 0+, 1" and 1+, 
respectively. For the (2, 2) system with 0^ symmetry, our 
results compare well with the value 2.5092 obtained by 
von Stecher and Greene [l| using the same approach as 
that employed here. For the (2,2) systems with 1~ and 
1+ symmetry, the Sq ^^j^ value has, to the best of our 
knowledge, not been calculated within the hyperspheri- 
cal coordinate approach. However, using scale invariance 
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FIG. 1: (Color online) Scaled hyperangular eigenvalue 
so,unit(i?) for the (2,2) system with /t = 1 at unitarity. Cir- 
cles show the scaled hyperangular eigenvalue so, unit (-R) as a 
function of ^/Jlro/R for (a) — 0"*" symmetry, (b) = 1^ 
symmetry, and (c) — l""" symmetry. The dotted lines show 
three-parameter fits (using a second-order polynomial). The 
solid horizontal line in panel (a) shows the result from Ref. [l|] 
for So^„it while the solid horizontal lines in panels (b) and (c) 
show the results from Ref. [31 for Sq u^jt (see the text for de- 
tails). The agreement between our results and those from the 
literature is very good. 



arguments [43|, Sg ^^j^ can be extracted from the energy 



cix + C2X , spectrum of the harmonically trapped (2, 2) system. The 



^0 unit values for the zero-range system at unitarity, ob- 
tained by analyzing the energy spectrum of the trapped 
4.5978 for = 1~ symmetry and 



system, are s'^^^^^ 



So.unit = 4.0820 for = 1+ symmetry Our values 
reported above are in very good agreement with these lit- 
erature values, indicating that the HECG approach is ca- 
pable of reliably describing strongly-correlated few-body 
systems with finite angular momentum and positive and 
negative parity. 
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FIG. 2: (Color online) Convergence of the scaled hyperangu- 
lar eigenvalue with increasing number of basis functions A^i,. 
Solid and dashed lines show the scaled hyperangular eigen- 
value So, unit (i?) as a function of A^;, for the (3, 1) system with 
1+ symmetry and k = 1 at unitarity for y^flro/R — 0.02 and 
y/JIro/R = 0.005, respectively. The inset shows a blowup. 



To illustrate the convergence of the HECG approach 
with the number of basis functions, we consider the (3, 1) 
system with l"*" symmetry (k = 1 and l/og = 0). Solid 
and dashed lines in Fig. [2] show the scaled hyperangu- 
lar eigenvalue so,unit(^) as a function of the number of 
basis functions Nb for ,/JIro/R — 0.02 and y/Jlro/R = 
0.005, respectively. For these calculations, we consid- 
ered Nt — 4800 and 6000 trial functions, respectively, 
for each new basis function selected. Figure [2] shows that 
the description of the system becomes more challenging 
as the separation of length scales increases, i.e., as the 
ratio y/JIro/R decreases. Moreover, it can be seen that 
So, unit (^) shows a few "shoulders", suggesting that there 
is room to improve upon the selection of the basis func- 
tions. Possible improvements may include gradient op- 
timization techniques, which have been successfully em- 
ployed in electronic structure calculations [il], or a re- 
fined trial and error procedure. Nevertheless, the HECG 
approach in its present implementation yields good con- 
vergence for basis sets consisting of around 100 basis 
functions for the systems under study. The so.unit(-R) for 
the largest basis set considered are 4.08194 and 4.0820 
for y/JIro/R = 0.02 and y^ro/R = 0.005, respectively. 
The corresponding basis set errors are estimated to be 
smaller than 0.00002 and 0.0002, respectively. The de- 
pendence of So. unit (^) on tq is relatively weak for the 
(3,1) system with 1"*" symmetry and we find the extrap- 
olated value s^i;^„ij = 4.0820(3). This result agrees well 
with the value of Sq unit = 4.0819 extracted from the trap 
energies [il ]. 

Circles in Fig. [3] shows the scaled hyperangular eigen- 
value So,unit(-R) for the (3,1) system with ^JJir^/R — 
0.005 and 1"*" symmetry (k = 1 and l/os = 0) as a func- 



FIG. 3: (Color online) Convergence of the scaled hyperan- 
gular eigenvalue with increasing number of grid points. The 
circles show the scaled hyperangular eigenvalue so,unit(i?), cal- 
culated for a basis set consisting of Nb = 95 basis functions, as 
a function of the order A'ordcr per sector used to perform the 
numerical integration over the angle 72; as discussed in the 
text, the numerical integration is divided into three sectors. 
For A'ordcr = 50 (not shown), the numerical integration breaks 
down (it yields a value that deviates by 10% from the exact 
value). The calculations are performed for the (3, 1) system 
with 1"*" symmetry, 1/as =0, k = 1 and ^JJlro/R — 0.005. 
The dotted line is shown as a guide to the eye. 



tion of the order A'ordcr used to perform the numerical 
integration over the hyperangle 72. As discussed above, 
the numerical integral is divided dynamically into three 
sectors, yielding a total of 3iVorder integration points. For 
the calculations shown in Fig. [3J we used a basis set with 
Nb = 95. Figure [3] shows that the scaled hyperangular 
eigenvalue so^unit(^) is, for the system considered, accu- 
rate to better than 0.1% for iVordcr ^ 60. It is important 
to note, though, that while iVordcr — 60 yields quite ac- 
curate results, iVordor = 50 yields completely unreliable 
results for the parameter combination and basis set con- 
sidered. In practice, we perform the optimization of the 
basis set for fixed iVorder- At the end of the construction 
of the basis set, we increase A'ordor to ensure that the 
results are independent of the integration scheme em- 
ployed. In general, the smaller the ratio ^/Jlro/R, the 
larger A'order (assuming the same L, 11, as and k). 

As pointed out in Sec. llVi the expressions for the ma- 
trix elements presented in the appendices apply to any 
number of particles. To explicitly confirm this, we per- 
formed calculations for the non-interacting five-body sys- 
tems with 0^, 1~ and 1+ symmetry and found hyperan- 
gular eigenvalues consistent with what is expected. Since 
the matrix elements for the two-body interactions are of 
the same form as those for the overlaps, treatment of the 
non-interacting systems suffices for testing the analytical 
expressions presented in this paper. The computational 
demands will, of course, increase as the interactions are 
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FIG. 4: (Color online) so,unit{R) as a function of ^fJlro/R 
for the (3, 1) system with 1+ symmetry at unitarity for k = 1 
(circles), k = A (pluses), k = 8 (squares), k = 9.5 (diamonds), 
and K = 10 (triangles). Dotted lines show three-parameter fits 
to the data. 

turned on. Assessing the performance of the outlined 
formalism for strongly-correlated five-body systems is be- 
yond the scope of this paper. 



VI. (3, 1) SYSTEM WITH 1+ SYMMETRY 

This section considers the (3, 1) system with 1+ sym- 
metry at unitarity for various mass ratios. The Sq y^jt 
value for these systems has been determined previ- 
ously m, 113] by investigating the (3,1) system under 
spherically symmetric harmonic confinement using the 
stochastic variational approach combined with geminal 
type basis functions, which are neither characterized by 
good angular momentum and corresponding projection 
quantum numbers nor good parity. As a result, the ear- 
lier calculations were restricted to comparatively large 
''o/flho values, where Oho denotes the harmonic oscillator 
length of the external confinement. The present work 
determines Sq ^j^j^ for various mass ratios k by employ- 
ing the HECG approach. For comparative purposes, 
we repeat the trap calculations using the standard EGG 
approach; however, instead of using geminal type basis 
functions we employ basis functions which are character- 
ized by good L, Ml and IT quantum numbers, thereby 
allowing us to reduce the basis set extrapolation error 
and to treat systems with smaller ro/aho than consid- 
ered earlier. 

Symbols in Fig. H] show so,unit(^) [see Eq. ([55])]. ob- 
tained by the HEGG approach, as a function of ^JJir^jR 
for K — 1,4,8,9.5 and 10. Dotted lines show three- 
parameter fits. The resulting Sq ^^^j^ values are summa- 
rized in column 2 of Table HI The errorbars are primarily 



TABLE I: The first column shows the mass ratio k. Columns 
2 and 3 show the Sojjnit values for the (3, 1) system with 
1^ symmetry at unitarity obtained by the HECG approach 
and from the extrapolated zero-range energies of the trapped 
system (see the text for details). For comparison, column 4 
shows the results from Ref.j47(]. The k = 1 entry in the third 
column is taken from Ref. [44!]. 
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due to the extrapolation to the y^ro/i? — > limit and 
only secondarily due to the basis set extrapolation error 
of the So, unit (^) for each R. Figure S] suggests that the 
leading order correction to Sg.unit proportional to 1/i? 
for the K values considered (k > 1). This is in agreement 
with what has been found analytically for the (2, 1) sys- 
tem with 1^ symmetry [4^. Figure |3] also shows that the 
range dependence of so,unit(-R) increases with increasing 
K. For larger k, the range dependence appears to be more 
complicated and the determination of the corresponding 
So, unit (^) values is beyond the scope of this paper. 

Gircles in Fig. [5] show the Sg unit values obtained from 
the HEGG approach as a function of k. For comparison, 
crosses show the Sg ^^j^ values obtained by extrapolat- 
ing the finite-range energies of the trapped system to the 
^'o/oho limit. These Sq^^j^ values are reported in 
column 3 of Table H] and are calculated following the pro- 
cedure discussed in Ref. [13] for equal masses. The agree- 
ment between the two sets of calculations is very good. 
Golumn 4 of Table U shows the Sg ^j^j^ values obtained 
earlier |43]; these earlier calculations were restricted to 
larger rg/aho and are less accurate than the values cal- 
culated in this work. 

Following the discussion of Refs. [H, [H, [H, 113, Elj , 
the Sq value indicates whether the system behaves 
universal or not. For two-component Fermi gases with 
zero-range s-wave interactions, e.g., the Sq ^^j^ value is 
larger than 1 and, correspondingly, the system proper- 
ties are fully determined by as- If Sg ^j^j^ < 1, the so- 
lution to the hyperradial Schrodinger equation, which 
is a second order differential equation, can — at least in 
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FIG. 5: (Color online) So,^nit a-s a function of k for the (3, 1) 
system with 1+ symmetry at unitarity. The circles show the 
■So.unit values determined using the HECG approach while the 
crosses show the Sq unit values determined by extrapolating 
the trap energies obtained by the standard EGG approach to 
the zero-range limit. The errorbars increase with increasing 
K. The agreement between the Sg u^jt values determined by 
the two different approaches is very good. The main panel 
and inset show the same data on different scales; the main 
panel and inset use respectively logarithmic and linear scales 
for K. 



principle — contain contributions of the "regular" and "ir- 
regular" solutions. If the irregular solution contributes, 
the system is said to behave non-universal since its prop- 
erties depend not only on the s-wave scattering length 
but additionally on a second parameter. Applying these 
arguments to the (3, 1) system with 1"*" symmetry and 
using that Sq unit > 1 mass ratios considered, the 

present study supports the finding that the four-body 
bound states found in Ref. j21|] for k > 9.5 and positive 
s-wave scattering length are universal. 



VII. CONCLUSIONS 

This paper extended the HECG approach, which had 
previously been formulated for and applied to three- and 
four-particle systems with = 0+ symmetry [H-ll], to 
states with 1~ and l"*" symmetry. The developed frame- 
work is applicable to systems with any n; realistically, 
though, applications in the not too distant future will 
likely be limited to systems with up to five particles. 
This paper emphasized a unified formulation for solving 
the hyperangular Schrodinger equation. In particular, 
many of the resulting equations apply to any particle 
number and symmetry, suggesting a numerical im- 
plementation in which most subroutines can be used for 
any particle number n and any symmetry; only a 
few subroutines specific to the values of n, L and 11 are 
needed. 



As a first application, we considered the (2, 2) and 
(3, 1) systems at unitarity. In particular, we solved the 
hyperangular Schrodinger equation for the energetically 
lowest-lying eigenvalue in the small ^/JItq / R regime and 
extracted the corresponding Sq unit values. Our results 
are in very good agreement with results from the liter- 
ature and with results determined by an alternative ap- 
proach. The Sq ^jjjj values for the (3, 1) system at unitar- 
ity with 1+ symmetry consisting of three heavy identical 
fcrmions and one light impurity particle are relevant to 
the (3, 1) system with positive a^. In particular, the re- 
sults obtained in this paper lend strong support that the 
bound states of the (3, 1) system with positive s-wave 
scattering length found in Ref. [2l| are universal, i.e., 
fully determined by a^. While we were able to reliably 
describe four-body systems for which the hyperradius R 
is 200 times larger than the scaled range yTflrQ of the two- 
body potential, pushing this ratio to much larger values 
may be challenging. 

In the future, it will be interesting to combine, as done 
in Refs. [IMS'], the developed framework with a standard 
R-matrix approach and to describe the scattering prop- 
erties of four-particle systems with finite angular momen- 
tum. The generalization of the developed formalism to 
states with other symmetries, which amounts to de- 
termining the corresponding d-, b-, p- and g-coefficients, 
is tedious but straightforward. Other possible extensions 
include the generalization of the approach to cold atom 
systems in a wave guide geometry. 



Appendix A: Definition of (symmetry-dependent) 
auxiliary quantities 

We first introduce a number of auxiliary quantities and 
then define quantities specific to the basis functions with 
= 0+, 1~ and 1"^ symmetry. 
The matrix S is defined as 



where 



and 



i£ — AbAb + AbAbi 



Ab - iUBfAUB 



A'b = [UBfA'UB- 

Similarly, we define 

B.B = Ab + Ab ■ 
We define the vectors uj^b and u'j b U — ^ ^^'^ 2), 



Uj.B = (ILb) 



and 



u'j^B = iU-Bfu'^ 



(Al) 

(A2) 
(A3) 
(A4) 

(A5) 
(A6) 
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Lastly, we define 

a{j,k) = ui^B{j)u[B{k), 

a{j,k) = U2,b(j>2 R(fc), 



(A7) 
(A8) 



Here, ^dsOj k) denotes the square of the matrix element 

ABij:k). 

The non-zero 6-coeflicients are 6'-'^^ , foj-^"* , 6^-^'' and b^^^"* , 



6(0) = -3Tr(B5), 



(A20) 



g{j,k,l) = a{j,j)a{l,k) + a{l, k)a{j, j) - 

a{l,j)a{j,k) ~ a{j,k)a{l,j), (A9) 



hU^k,l) ^g{l,j,k) +g{l,j,k), 



(AlO) 



(All) 



k, I, m) — a{j, k)a{l, m) + a{k, j)a{l, m) + 

a{i, k)a{m, I) + a{k, j)a{m, I), (A12) 

G(j, k, I, m) = a{i, l)a{m, k) + a(j, ■m)a{l, k) + 

aik, l)a{'m, j) + a{k, 171)0,(1, j), (A13) 



and 



f{j,k,l,m) = F{j,k,l,m) + F{l,m,j,k) - 

G{j,k,l,m)^G{l,m,j,k). (A14) 



Sections lA lllA 31 give explicit expressions for the d- 
, p-, q- and 6-coefficients for the basis functions with 
0"*", 1~ and 1+ symmetry, respectively. In what follows, 
the elements of the vector Uj^s are denoted by Uj_B{k), 
k = I,-- - ,11 — 1; the same notation is adopted for the 
elements of other vectors. The elements of the matrix 
S_ are denoted by S_(k,l) with k and I — 1, • • • ,n — 1; 
the same notation is adopted for the elements of other 
matrices. 



1. 0"*" symmetry 



The only non-zero d-coefficient is d^^\ 
= 1. 

(2) 

The only non-zero p-coefficient is , 

.(2) 



The non-zero (/-coefficients are ^j^'' and 

(2) n-i (2) 



(A15) 
(A16) 

J22) 



(A17) 



and 



(22) 



^ABij.MBik,k) + -Alij,k) 



(A18) 



i?"^(Al9) 



= /?! ~ m j) + (3n - 3)/3,i?-2, (A21) 



= + 2Ab ij, Mb j)] (A22) 



and 



1,(22) 



-2p,l3uR-'' + 2As(j, jU^(fc, k)R-^ + 



2As{k, k)A:sU,j)R-' + -AsU, k)A!sij, fc)i?-2.(A23) 



2. 1 symmetry 

The only non-zero d-coefhcient is d^^) , 

df=aU,j). (A24) 
The non-zero p-coefficients are p^^^ , p^f^ and p^^? , 

pf =-a{j,3)R-\ (A25) 



and 



P^-^ ^Ab{j.3)<J.3)R-\ 



AB{k,k)a{j,j)R-^ + 



(A26) 



(A27) 



The non-zero g-coefhcients are q] ' , 9] , q) , q 



(4) ^(22) ^(6) (24) 



and grj^p. 



qf=5ir^pf, 



(22) oD-1 (22) 

ij,k -^R Pj,k^ 



if = -Al{j,j)a{j,j)R- 



(A28) 



(A29) 



(A30) 



(24) 



-A|(fc,fc)a0-,j)i?-2- 

^A|(j,fc)a(fc,fc)i?-2- 
^Asik, k)As{j, k) [a(j, fc) + a(fc, j)] i?-^ (^31) 
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and 



3. l"*" symmetry 



„(222) _ 



[aij,k) + a{k,j)] RXki2) 



The non-zero 6-coefficients are 6^-^' , b'^^^ , b'^p , b^^"^^ , b'^^^ 



and b 



(222) 



The only non-zero d-coefficient is d^f'^\ 

4^fe^ = "■{3d)o.{k, k) + a{k, k)a{j,j) - 

a{j, k)a{k,j) - a{k,j)d{j, k). (A39) 



The non-zero p-coefficients are pj^^'' , J'j^fe and pf^'^i , 



„(22) _ ,p-lj(22) 

Pj,k - "-j.k ' 



(A40) 



bf = -2{3n-4)a{j,j)R-^- 

n-l 

2j2[ABij,k)a{k,j)+A's{.hk)a{j,k)] - 

3Tr{Bs)a{j,j), (A33) 



fe=i 



bf = - S{j,j) + (3n - l)l3jR-'] a{j,j), (A34) 



and 



P^^'=R-'d^UBik,k), 



pf,S = R~'€kAB{i,i) + 

lR-'ABU,k)h{j,k,l). 



-lj(22) 



(A41) 



(A42) 



= + '^AbU^MbUJ)] a{j,3)R-\ (A35) 



The non-zero g-coefhcients are qf^ , Qj^k^, ifkf^ 1j\k^' 



(44) (224) , (2222) 
9j,fe ' 1j,k,l l3,k,l,m' 



[Pj - Sij,j} + (3n - l)f3jR-^] a{k, k) + 
[Pi - S{k, k) + (3n - aO; j) - 

\ m, k) + S{k,j)] [a{j, k) + a{k,3)\ , (A36) 



(22) D_i (22) 

ij,k-R Ak' 



„(24) _ 5D-1„(24) 
1j,k — "^-^ Pj,k ' 



(A43) 



(A44) 



b?t^ = 



[-01 + '^Aeik, fcU'sCfc, k)] a{j,j)R-^ + 
[-2/3,/3fc + 2As(i,j)A^(fc,fc) + 

2Aj,{k, k)A's{j,j) + ^AbU, kJA'sU, k)] x 
a{k,k)R-^ + 

^ [AbU, k)A'B{k, k) + AB{k, kUeij, k)] x 



[a(j, k)+a{k,j)] R~ 



(A37) 



„(222) _ rD-i„>,^^^; 

^j,k,l — Pj,k,l ' 



1J222) 



and 



<ifP=-df^Al{k,k)R-\ 



& = -2df^HBiJ,j)ABik,k)R-' - 
ld^H%ij,k)R-^ 



(A45) 



(A46) 



(A47) 



,(222) 



[-2(3j(3k + 2As{j,j)A!s{k,k) + 

2AB{k, k)A'B{jJ) + IAbU, k)A'{3, k)] x 
a{l,l)R-'^ + 

^ [AB{j,k)A'B{lJ)+ABil,l)^BU,k)] X 

[a{j,k) + aik,j)] R-^ + 

I [ABij,l)A:B{k,l)+Asik,l)A!sij,l)] X 



[a{j,k)+a{k,j)]R- 



(A38) 



„(224) _ 
'^j,kj ~ 



AB{j,j)AB{l,l) + ^AUj,l) 



2R-^S^f - 



AB{k,k)Aj,{l,l) + ^Al{k,l) 



2R-H^? - 



A%{l,l)R-^d 



'j,k 



^[ABij,k)As{l,l)h{j,k,l) + 

ABU,i)AB{iJ)KiJ,k) + 

AB{k,l)AB{l,lMk,l,3)]R-\ 



(A48) 
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and 



2As(/,/)^B(m,m) + -A|(/,to) 

2ABij,j)ABik,k) + ^Alij,k) 



(2222) ^ 
^j.k,l.'m 



3 [AbU^ 'm)h{j, k, I) + 

AB{3,k)ABil,lMj,k,m)]R-^ - 
4 

- [Asil, m)As{k, k)h{l, m,j) + 
AB{l,m)AB{j,j)h{l,m,k)]R-'^ - 

- [Asil, k)As{k, m)h{l, m,j) + 

Aaih j)AB{h'^W^'rn,k)\R-'^ - 
-^[ABi3,'ni)Ag{m,k)h{i,k,l) + 
ABU,l)ABii,k)h{j,k,m)]R-^ - 
I IAbU, k)ABil, m)] f{j, k, I, m)R-\AA9) 



The non-zero 6-coefficients are 



ri-l 

2 J2 [^bU, l)9{k,j, I) + A!s{j, l)g{k,j, + 

Aeik, l)g{j, k, I) + A^B (fc, l)g{j, k, I) - 
A4f]' (A50) 



(22) 



-2 ,(22) 



= [Pi - S{k, k)] d}f^ + /3,(3n + l)i?-24;'^^(A51) 



^(22) 



2j(22)/ 



€? = [A^-^(M)]4?- 

\h{3, k, l)S{j, k) + A(3n + l)i?-'4?' (^52) 



blf = [ - Pjf^k + ABi3,j)A'Bik, k) + 
4b (fc, A)4b(j, j) + \Ab{3, k)A's{j, k)\2R-Hf^\kbi) 



C = [-f^l + 24B(fc, k)A!^{k, k)] R-^dfl\ (A54) 



^AB{3MB{3,l)\2R-%y + 
[ - Ml + Aeik, k)A:B{l, + Ab{1, 04b (fc, k) + 

^ABik,l)^Bik,l)]2R-^dfP + 
[^Pf + 2As{l,l)A!sil,l)] R-^df;t^- 
\ [Ab{3, k)AB{l, I) + 4b (j, fcUB(«, I)] R-^h{j, k, I) - 

I [AbU, i)Ab{i, + A'bU, i)A'b{i, I)] R-^Kh h k) - 

I [ABik, l)As{l, I) + A'eik, l)A'sil, I)] R-^h{k, 1,^5) 



and 



,(2222) ^ 

-{2AB{l,l)AB{m,m) + 2A'B{l,l)A!Bim,m) + 
l[Al{l,m)+A:l{l,m)] }ii-^4f - 
{24b(j, jUB(fc, k) + 2A!sU,j)A:B{k, k) + 
l[Al{3,k)+A:i{3,k)] }R-'d^^^- 

4 

3 [Ab U, k)AB (m, m) + 4s (i> ^)4s {m,m)] x 

h{j,k,l)R-^- 

\ [Ab{3, k)AB{l, I) + ^b{3, k)A'B{h 0] X 
h{j, k,m)R~^ — 

^ [ABil,rn)AB{k,k) + A!B{l,rn)A!B{k,k)] x 

h{l,m,j)R-^- 

^ [4s(^."^)4B(i,i)+4BG.m)4B(i,j)] x 

h{l,m,k)R-^ - 
- [AB{l,k)AB{k,m) + A:B{l,k)A;s{k,m)] x 

h{l,m,j)R-^- 

g [4b {1,3)Ab U, m) + Kb {1,3)A!b (j, m)] x 

h{l,m,k)R-^ - 
g [4B(i>TO)4BK^)+4s(j,m)4^(m,fc)] X 

h{j,k,l)R-^- 

[4b (i, 04b k) + A:b{3, 04b k)] h{j, k, m)R-^ - 

g 

g [4B(i.^)4B(^"i)+4B(i,fc)4B(^"^)] X 

f{j,k,l,m)R{AL%) 
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Appendix B: Sketch of derivation of matrix elements 
for n — 3 and = 1^ symmetry 



This appendix derives the overlap matrix element for 
n = 3 and = 1~. To evaluate the overlap matrix 
element, we write, using Eq. (PU]) . 



ip{A', u[, x)'il;{A, ui,x) = 3v[ 3W1.3 exp 



Bx 



,(B1) 



where B_ is defined in Eq. (j36p . Using the transformation 
y — U_gX, we find 

a(l, 1)2/2^3 + [a(l, 2) + a(2, 1)] yi^3y2,3 + a(2, 2)^1, 3,(B2) 



where a{i,j) is defined in Eq. (|A7P and where j/j^3 denotes 
the z-component of the vector yj. Using j/i 3 — yi cosi?! 
and j/2,3 = 2/2 cos in Eqs. (|B1I) and (IB2p . we have 



V-V = 3{a(l,l)yJ cos^i?! + 
[a(l, 2) + a(2, 1)] yij/a cost^i cosiJz + a(2, 2)y^ cos^ i^z} x 



exp 



-^{f3iyl + f32yi)lm 



The derivation sketched above for the overlap matrix 
element can be fairly straightforwardly generalized to ar- 
bitrary n. To calculate the matrix element {'ijj\TQ\tp)\H, 
we use Eq. (|20p . i.e., we separately calculate the matrix 
elements {ip\T,-ci\'ip)\R and {ilj\Tii\ip)\ji. The evaluation of 
these matrix elements is not fundamentally difficult but 
somewhat tedious and lengthy. 



Appendix C: Definitions of sc-coefiicients, Afi and M2 



The sc-coefficients entering into fo are given by 



scoo 



where the /3j denote, as before, the eigenvalues of the 
matrix B. When integrating over yij/2, the cross term 
averages to zero and we have 



i=3 

n-l 

/ . ["'j,k VjVk + "-j^k VjHk 
k>j;j,kyil.,2 
n-l 

/ , ["■j,k yjVk^o,j,k VjVk 

j=3,k=3,k^j 
n-l 

E [ 

k>r,l^j,k;j,k,lyil,2 
n-l 



,(222) 2 2 2 I ,(224) 2 2 4 

dj.k/VjVkyi + d)^k,iyjykyi 



E 



dTkZy'ylkm) 



ky j :ly j :myl;j^k^l^m-j ,k,l ,m^l ,2 



2tt rl 1-2-K i-l 



-IJO J-1 



(ip'ip) \ Rd COS iSidipid COS '&2d^2 = 
(47r)2 [a(l,l)y2 + a(2,2)y2] x 



exp 



-^(/3iy? + /32yi) 



(B4) 



(2) 



Comparison of Eq. (|B4[) with Eq. shows that d' 
a{j,j) and that all other ^-coefficients are zero. This 
agrees with the expressions given in Appendix lA 21 

Next, we consider the integration over 71. To this end, 
we replace yi and y2 in Eq. (jB4p by i? sin 71 and i? cos 71 , 
respectively. According to the discussion in Sec. IIVI 71 
can take values between and 7r/2. Multiplying both 
sides of Eq. (jB4p by sin^ 71 cos^ 71 [see Eq. ([55]) ] and in- 
tegrating over 71, we find 



(V'V)|fld^f1= (47r)2_exp 



{[a(l,l) + a(2,2)]i?2/i(C) + 
ha(l,l) + a(2,2)]i?2/2(^)}(B5) 

Applying the definitions from Appendix [X] and [Cl it can 
be verified that Eq. (jBSP agrees with Eq. (|46l) [note that 
Eq. does not contain a factor of (47r)^ while Eq. (jB5[) 
does; the reason is that fo is defined without this prefac- 
tor]. 



SC20 



n— 1 



3=3 

n-l 

V Tr/^^^^^y^ 2 ^(224) 2 4 

j,k=3;k^j 

n-l 

/ . ^j,k,iyjyk + 

j<k-j,k^l,2 
n-l 

E '^ij.k.jy'jyki^f' 

l>k;i^k^l;i ,k,l^l,2 



n— I 



3=3 



n-l 



/ . ["2j,fe + «2j,fe yjj^fe 

j,fe=3;/c5iij 
n-l 

/ . "■j,k,2yjyk + 

j<k-j,k^l,2 
n-l 

/ . "'2,j,k,iyjykyi ^ 

l>k;j^kyil;j,k,l^l,2 



(C2) 



(C3) 
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n-l 



SC40 



(24) 2 



+ 



i=3 



2^ "'j,k,iyjyk, 



i<fc;j,fc#l,2 



(C4) 



SC26 = 7 



(C12) 



and 



SC62 = (C13) 



SCoi 



(4) 



n-l 



+ 



3=3 



,(6) , ,(26) 2 

J=3 



n-l 

J=3 



SC22 



^(22) 
"1.2 



(C5) 



(C6) 



(C7) 



n — 1 

Er r ,(222) , ,(222) , ,(222)1 2 , ,(224) i\ , 
I [di,2j + rfi j,2 + ^^2j,i J yj + di,2jyj \ + 



3=3 

n-l 

E,(2222) 2 2 

i<A;u\fe/1.2 
n-l 

/ . "ij2fe Vjyk^ 

j#fcU,Ml,2 



(C8) 



Equation (|Cip - (IC13l) also apply to /p, fq and /o if the 
d-coefRcients are replaced by the p-, q- and 6-coefRcients, 
respectively. The quantities Mi and M2 are defined 
through 



Ml 



[2SC00 + (SC20 + SC(32)HR^ + 



(sc4o + sco4){HR^f + (sceo + sco6)(i/i?')'] C (CM) 



and 



M2 = C"^Maux + 4 [ - 2SC40 - 2SC04 + 2sC22 + 

(-3SC60 - 3SC06 + SC24 + SCA2)HR^ + 
(sC26+SC62)(i/i?')'](ffi?')' + 

[ — SC20 + SC02 + (— SC40 + scoi)HR^ + 

(-SC60 + sco6)(i/i?')']i?i?\,(C15) 



where 
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Maux = -y(-SC26 + SC44 - SCG2)iHRYC^ + 

3 [ — SC60 + SC06 — SC24 + SC42 + 

(-SC26 + SC62)i?i?'] (HR^f. {C16) 



SC4:4: 12 ' 



,(24) , ,(224) 2 

SC24 = di,2 + 2^aij,2 > 

J=3 



n-l 

,(24) , ,(224) 2 

SC42 = d2,l + Z^«2j,l 
J=3 



(C9) 



(CIO) 



(Cll) 
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